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We performed three-dimensional Particle-in-Cell simulations of magnetic reconnection with mul¬ 
tiple magnetic null points. Magnetic field energy conversion into kinetic energy was about five times 
higher than in traditional Harris sheet configuration. More than 85% of initial magnetic field energy 
was transferred to particle energy during 25 reversed ion cyclofrequencies. Magnetic reconnection 
in the cluster of null points evolved in three phases. During the first phase, ion beams were ex¬ 
cited, that then gave part of their energy back to magnetic field in the second phase. In the third 
phase, magnetic reconnection occurs in many small patches around the current channels formed 
along the stripes of low magnetic field. Magnetic reconnection in null points presents essentially 
three-dimensional features, with no two dimensional symmetries or current sheets. 


Null-point magnetic reconnection is thought to be the 
main source of energy release in solar flares, Earth mag¬ 
netosphere and other astrophysical plasmas. Simulations 
of 3D null point reconnection in magnetohydro dynamic 
(MHD) approach were reported in [IHS]; null-point re¬ 
connection regimes were reviewed and classified by [4]; 
in the papers 3D null point reconnection was stud¬ 
ied in kinematic approach; [ZHg investigated particle ac¬ 
celeration in the vicinity of 3D null-points. In many as¬ 
trophysical applications, such as planet magnetospheres 
or solar flares, fluid models are not sufficient to describe 
magnetic reconnection unmi]. To our knowledge, only 
m reported Particle-in-Cell (PIC) simulations of an iso¬ 
lated null point in solar corona. To understand the basic 
kinetic and 3D effects of null point reconnection, we have 
performed PIC simulation of a cluster of null points in 
the simplest possible cofiguration that excludes the com¬ 
plications due to boundary condition or external sources. 
We have found that energy release in such system is much 
more efficient than in traditional Harris sheet configu¬ 
rations. Three-dimensional effects dominate in the dy¬ 
namics throughout the simulation, and no current sheets 
or other 2D structures are formed during reconnection. 
We believe that both laminar and turbulent reconnection 
regimes were realised in our experiment, characterized by 
different energy dissipation and particle acceleration [13- 

[ig. 

Making the very first step towards kinetic models of 
3D null point reconnection, we designed our simulations 
to be as simple and straightforward, as possible. In par¬ 
ticular, we impose boundaries in our simulation domain 
periodic in all directions. In a fully periodic box, it can 
be proved that the minimum number of points where 
magnetic field vanishes, is 8. We have chosen a non¬ 
equilibrium initial state with uniformly spaced 8 null 
points: 
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where 5o is the magnetic field amplitude; Ly, and Lz 
are the sizes of the simulation domain in the correspond¬ 
ing directions. It is easy to show that this configuration 
satisfies the condition V • B = 0. Besides 8 null points, 
the condition 5 = 0 also holds along 9 “null lines” lied 
up in the planes Y = irn, where n is an integer (Fig. [^. 

The collisionless plasma simulations were carried out 
using iPIC3D, fully kinetic electromagnetic PIC code 
with implicit time stepping m- Simulation domain rep¬ 
resents a cubic box with dimensions L^ ^ Ly x Lz = 
IDdi X ID di X IDdi with 256^ cells and 64 particles of 
each specie per cell. Our plasma consisted of electrons 
and ions with mass ratio mijm^ = 25 and tempera¬ 
ture ratio TijT^ = 5; the particles were initialized with 
Maxwellian velocity distribution with electron thermal 
velocity Vthejc = 0.0346. The initial particle density was 
uniform with no = I; the initial magnetic field amplitude 
was 5o = 0.03. It is important to note, that with such 
parameters the electron spatial scales are well resolved 
in our simulation: the electron skin depth de = 0.2 di = 
5Ax, where Ax is the grid spacing. The time step was 
set relative to ion plasma frequency: At = 0.125/(1;^^. 
For consistency, in the following all times are normal¬ 
ized to ion cyclofrequency computed with the initial field 
amplitude Vtci = eBoIrriiC = 0.052 

Our simulation domain initially includes null points of 
both A and B types m, as shown in Figure All the 
null points are interconnected with each other by sepa¬ 
rators, formed by the intersections of separatrix surfaces 
and spines. Many of these separators are the artifacts of 
the highly symmetric configuration, and must not neces¬ 
sarily exist in nature [18] ; they are destroyed very quickly 
after beginning of the simulation. 
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FIG. 1. a, b, c: Snapshots of magnetic field at different times (labelled in the upper left). Field lines are plotted in different 
colors; each color corresponds to the lines passing by one of 8 initial null points. Isocontours of low magnetic field B = O.IBq 
are plotted in grey, d: Snapshot of current J made at Qcit =1.6. 


At Qcit = 1-6 (panel b in Fig. no aforementioned 
“null line” is visible, and the connectivity of many field 
lines has changed. Initially, gas and magnetic pressures 
are not balanced, and the magnetic fields surrounding the 
“null lines” are compressed by plasma. Particles move 
towards these regions and establish the current channels 
that dominate the evolution of the system, as shown in 
panel d of Figure The energy released by magnetic 
field in this phase of the simulation is transferred mainly 
to ions accelerated in the current channels. The radius 
of the cross-section of such channel is of the order of di. 
Current channels become unstable later in time, they are 
continuously distorted, and multiple stripes of small B 
intensity are formed along them. As we will demonstrate 
below, the reconnection sites are associated with these re¬ 
gions. There is no initial current through the null points 
in our configuration, and the null points do not move 
during the simulation. They are destroyed in pairs (see 


article supplemented material). 

At the end of the simulations all 8 null points are de¬ 
stroyed, and the magnetic field organization is not rem¬ 
iniscent of the initial configuration. The magnetic field 
lines spiral around current channels established in the re¬ 
gions of low B. Note, that we have terminated the simu¬ 
lation at time = 26 when most, but not all the mag¬ 
netic field lines are reconnected. All the features forming 
during the simulation are essentially three-dimensional, 
and current sheets or other 2D structures are not created. 

Figure shows the evolution of the magnetic field en¬ 
ergy Wmagi total kinetic energy of the particles Wkim 
bulk and thermal kinetic energy for species {Wi)uik,ei 
Wbuik,i^ Wtherm,e, and Wtherm^i, Correspondingly). The 
electric field energy is two orders of magnitude smaller 
than other components of energy, and is not shown in the 
plot. Initially, the magnetic energy contributes 85% to 
the total energy of the system; particles contribute only 
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FIG. 2. Evolution of different components of energy (total 
kinetic energy of particles, magnetic field energy, bulk and 
thermal energy of both species). 

12%. At the end of the simulation, particle energy con¬ 
stitutes 87% of the total energy: 86% of magnetic energy 
is released within 26Q~-^\ This corresponds to decay rate 
A = {AWmag/^t) jWmag ~ O.OSSflci, which is almost 5 
times higher than A 0.0075f2ci observed in “classical” 
two-dimensional Harris current sheet reconnection setup 

[niEo]. 

In the beginning of the simulation, when the current 
channels are established, most of the released magnetic 
energy goes to the energy of ion bulk motions W^uik,!- 
Very soon, though, these motions are slowed down, and in 
the later phase of reconnection the temperature of species 
is increased rather than fluid motions. At the end of the 
simulations temperature ratio reduces to T^/Tg ~ 3, and 
fluid motions constitute only a few percent to the total 
energy in the system. 

As in mnn, we compute the electron-frame dissipa¬ 
tion measure 

Dq = J • (E — [Jg X B] /pe) + PgUg • E, (1) 

where J and Jg are the total and electron current densi¬ 
ties, correspondingly, and Ug are electron density and 
speed, and pc = e{ni — Ug) is charge density. The dissipa¬ 
tion De characterizes the work done by fields on particles 
J-E with corrections for convective motions J-[JgXB]/pg, 
and for charge separation pgUg • E. The latter quantity 
is by orders of magnitude smaller than two other compo¬ 
nents of De in our simulations. The first two components 
balance each other so that De is very small. 

For additional illustration, we have selected the 
“lower” half of our computational domain that spans 
from Z = 0 to Z = 5 di^ and contains 4 null points and 
sections of 6 “null lines”. Figureshows snapshots of De 
and parallel electric field = E • B/H in this slab. 

The magnetic reconnection dynamics may be divided 
into three distinct phases according to the energy evolu¬ 
tion. The first phase of rapid acceleration of particles and 
very high dissipation De lasts till Og^t ~ 2. During this 
time powerful currents are established in the low B value 
channels. The speeds of particle beams in these channels 


are up to 3.3Va for electrons and 1.6Va for ions, where 
Va is Alfven speed. Most of the electromagnetic energy 
at this stage goes to accelerate ion fluid motions; electron 
acceleration is much smaller. The dissipation regions and 
the regions of enhanced surround the current channels 
(left column in Fig. [^. 

During the second phase of magnetic reconnection, 
2.5 ^ Qeit ^ 7, work of the electromagnetic field on 
particles is negative, i.e., particles give some of their 
energy back to the field. The main contributor of this 
“reverse” energy exchange is the deceleration of ion cur¬ 
rents. Dissipation areas surround the current channels, 
while the regions of large values of E\\ are scattered ran¬ 
domly throughout the domain, so there is no correlation 
between these locations (middle column in Fig. H). The 
value of T^/Tg ^ 3 establishes in this phase, and does not 
substantially vary later. 

In the last phase of magnetic reconnection, small-scale 
reconnection events are chaotically distributed in the re¬ 
gions of low B (right column in Fig. [^. The locations of 
enhanced De and E\\ show well correspondence to each 
other. This phase is characterized by a steady growth 
of particle energy and a decrease of the magnetic energy 
(Fig.§. The released magnetic energy accelerates par¬ 
ticle random motions and increase their effective tem¬ 
perature. The steady (turbulent) reconnection process 
continues until all the null points are destroyed and all 
field lines reconnect (not shown in this paper). 

To summarize, we have performed PIC simulations of 
magnetic reconnection in a 3D domain that contains 8 
null points. The simulation was terminated at Qeit = 26, 
when all the initial null points vanish, and the major¬ 
ity of magnetic field lines spiral around the regions of 
low magnetic field and powerful electric currents. Dur¬ 
ing this period, 86% of initial magnetic field energy was 
transferred to particle kinetic energy, and the decay rate 
was almost 5 times higher than in two-dimensional Har¬ 
ris sheet reconnection. All the magnetic and current fea¬ 
tures that formed during the simulations were essentially 
three-dimensional, no current sheets or other 2D struc¬ 
ture were identified. 

Three phases of magnetic reconnection could be dis¬ 
tinguished: (i) a rapid establishment of powerful current 
channels, when most of energy is gained by ion beams; (ii) 
a negative dissipation phase, when ion beams slow down 
and return some energy back to electromagnetic field; 
(iii) a “steady” reconnection phase, characterized by mul¬ 
tiple small-scale reconnection events localized around the 
current channels. The first two phases lasted for ftei^ ~ 2, 
and Cleit ~ 5, correspondingly, while the third one was 
still developing when the simulation was terminated. Fi¬ 
nally, most of the magnetic energy was converted to the 
thermal energy of the particles, so that ion/electron tem¬ 
perature ratio changed from 5 in the beginning to 3 in 
the end of simulations. 

Interestingly, during the second phase the topology of 
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FIG. 3. Snapshots of the electron-frame dissipation measure De (upper panels) and parallel electric field E\\ (lower panels) in 
the lower half of the computational domain spanning from Z = 0 to Z = 5 at different times (labelled in the upper left of 
each column). Isocontours of low magnetic field B = 0.1 Bo are plotted in grey. 


magnetic field was changing, but the dissipation regions 
did not coincide with the regions of enhanced parallel 
electric field. During the “steady” reconnection, how¬ 
ever, these locations showed good correspondence to each 
other, surrounding the regions where the reconnection 
took place. 
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